home *** CD-ROM | disk | FTP | other *** search
/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc A) / Acorn User China CD-ROM (UK) (Disc A).bin / HENSA / MATHS / JACAL.ARC / jacal_ext < prev    next >
Encoding:
Text File  |  1993-10-03  |  6.5 KB  |  205 lines
  1. ;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
  2. ;;; Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
  3. ;;; See the file "COPYING" for terms applying to this program.
  4.  
  5. ;;; An algebraic extension is the root of a polynomial with more than
  6. ;;; one distinct value.  These values are not linked;  the difference
  7. ;;; between two algebraic extensions which are roots of identical
  8. ;;; polynomials is not 0.  Radicals have an additional rule that
  9. ;;; exponents of "positive" radicands commute.  For instance:
  10. ;;; (x^2)^(1/2) ==> x.  Notice that ((-x)^2)^(1/2) ==> x also.
  11. ;;; (-x^2)^(1/2) ==> (-1)^(1/2)*x.
  12.  
  13. ;;; algebraic extensions
  14. ;;; we want to find all extensions used by this poly except this poly.
  15. (define (alg:exts poly)
  16.   (let ((elts '()))
  17.     (poly:for-each-var 
  18.      (lambda (v)
  19.        (let ((er (extrule v)))
  20.      (if (and er (not (eq? er poly)))
  21.          (set! elts (adjoin v elts)))))
  22.      poly)
  23.     elts))
  24.  
  25. (define (application? v)
  26.   (and (not (extrule v)) (pair? (var:sexp v))
  27.        (not (eq? 'differential (car (var:sexp v))))))
  28.  
  29. ;;; we want to find all functionals used by this poly except.
  30. (define (var:funcs poly)
  31.   (let ((elts '()))
  32.     (poly:for-each-var 
  33.      (lambda (v)
  34.        (if (application? v)
  35.        (set! elts (adjoin v elts))))
  36.      poly)
  37.     elts))
  38.  
  39. ;;; algebraic and applications
  40. (define (chainables poly)
  41.   (let ((elts '()))
  42.     (poly:for-each-var 
  43.      (lambda (v)
  44.        (let ((er (extrule v)))
  45.      (if (or (and er (not (eq? er poly))) (application? v))
  46.          (set! elts (adjoin v elts)))))
  47.      poly)
  48.     elts))
  49.  
  50. ;;;alg:vars returns a list of all terminal vars used in this or in extensions
  51. ;;;used in this.
  52. (define (alg:vars poly)
  53.   (let ((deps '()))
  54.     (poly:for-each-var
  55.      (lambda (v)
  56.        (if (and (not (extrule v)) (null? (var:depends v)))
  57.        (set! deps (adjoin v deps)))
  58.        (set! deps (union (var:depends v) deps)))
  59.      poly)
  60.     deps))
  61.  
  62. (define (alg:square-free-var p var)
  63.   (alg:/ p (alg:gcd p (alg:diff p var))))
  64.  
  65. ;;; This is for equations
  66. ;;; Don't simplify a rule with itself
  67. (define (alg:simplify p)
  68.   (let ((exrls (map extrule (sort (alg:exts p) var:>))))
  69.     (if (memv p exrls)
  70.     p
  71.     (reduce-init poly:prem p exrls))))
  72.  
  73. (define (alg:clear-denoms p)
  74.   (do ((v (poly:find-var-if? (rat:denom p) potent-extrule)
  75.       (poly:find-var-if? (rat:denom p) potent-extrule))
  76.        (oldv "foo" (car v)))
  77.       ((not v) p)
  78.       (if (eq? (car v) oldv)
  79.       (eval-error 'could-not-clear-denominator-of:- p))
  80.       (set! p (alg:simplify
  81.            (poly:* p (alg:conjugate (rat:denom p) v))))))
  82.  
  83. ;;; This generates conjugates for any algebraic by a wonderful theorem of mine.
  84. ;;; 4/30/90 jaffer
  85. (define (alg:conjugate poly extpoly)
  86.   (let* ((var (car extpoly))
  87.      (pdiv (univ:pdiv extpoly (promote var poly)))
  88.      (pquo (car pdiv))
  89.      (prem (cadr pdiv)))
  90.     (if (zero? (univ:degree prem var))
  91.     pquo
  92.       (poly:* pquo (alg:conjugate prem extpoly)))))
  93.  
  94. ;;;This currently works only for univ extpoly
  95. (define (alg:mod poly extpoly)
  96.   (let ((p (poly:prem poly extpoly)))
  97.     (if (and (rat? p) (pair? extpoly)
  98.          (pair? (rat:denom p)) (eq? (car extpoly) (car (rat:denom p))))
  99.     (poly:prem
  100.      (poly:* p (alg:conjugate (rat:denom p) extpoly))
  101.      extpoly)
  102.     p)))
  103.  
  104. ;;; This section attempts to implement an incremental version of
  105. ;;; Caviness, B.F., Fateman, R.:
  106. ;;; Simplification of Radical Expressions.
  107. ;;; SYMSAC 1976, 329-338
  108. ;;; as described in
  109. ;;; Buchberger, B., Collins, G.E., Loos, R.:
  110. ;;; Computer Algebra, Symbolic and Algebraic Computation. Second Edition
  111. ;;; Springer-Verlag/Wein 1983, 20-22
  112. ;;; This algorithm for canonical simplification of UNNESTED radical expressions
  113. ;;; also has the convention that (s * t)^r = s^r * t^r.
  114. ;;; If the variable LINKRADICALS is #f then a new multiple value expression
  115. ;;; is returned for each radical.
  116.  
  117. ;;; this is actually alg:depth
  118. ;(define (rad:depth imp)
  119. ;  (let ((exts (alg:exts imp)))
  120. ;    (if (null? exts)
  121. ;    0
  122. ;      (+ 1 (apply max (map (lambda (x) (rad:depth (extrule x))) exts))))))
  123.  
  124. ;;; Integer power of EXPR
  125. (define (ipow a pow)
  126.   (if (not (integer? pow)) (math:error 'non-integer-power?- pow))
  127.   (cond ((expl? a) (if (< pow 0)
  128.                (make-rat 1 (poly:^ a (- pow)))
  129.              (poly:^ a pow)))
  130.     ((rat? a) (if (< pow 0)
  131.               (make-rat (ipow (rat:denom a) (- pow))
  132.                 (ipow (rat:num a) (- pow)))
  133.             (make-rat (ipow (rat:num a) pow)
  134.                   (ipow (rat:denom a) pow))))
  135.     (else (if (< pow 0)
  136.           (app* (list $ 1 (univ:monomial -1 (- pow) $1)) a)
  137.         (app* (univ:monomial 1 pow $1) a)))))
  138.  
  139. (define (^ a pow)
  140.   (cond
  141.    ((not (rat:number? pow)) (deferop _^ a pow))
  142.    ((eqn? a) (math:error 'Expt-of-equation?:- a))
  143.    (else
  144.     (set! pow (expr:normalize pow))
  145.     (let ((tmp #f)
  146.       (expnum (num pow))
  147.       (expdenom (denom pow)))
  148.       (cond
  149.        ((eqv? 1 expdenom) (ipow a expnum))
  150.        (linkradicals
  151.     (set! a (expr:normalize a))
  152.     (cond ((expl? a) (ipow (make-radical-ext a expdenom) expnum))
  153.           ((not (rat? a)) (math:error 'Non-rational-radicand:- a))
  154.           ((rat:unit-denom? a)
  155.            (ipow (make-radical-ext (poly:* (denom a) (num a)) expdenom)
  156.              expnum))
  157.           (else (ipow (make-rat (make-radical-ext (rat:num a) expdenom)
  158.                     (make-radical-ext (rat:denom a) expdenom))
  159.               expnum))))
  160.        (else
  161.     (app* (cond ((> expnum 0)
  162.              (set! tmp (univ:monomial -1 expdenom $))
  163.              (set-car! (cdr tmp) (univ:monomial 1 expnum $1))
  164.              tmp)
  165.             (else
  166.              (set! tmp (univ:monomial
  167.                 (univ:monomial -1 (- expnum) $1)
  168.                 expdenom
  169.                 $))
  170.              (set-car! (cdr tmp) 1)
  171.              tmp))
  172.           a)))))))
  173.  
  174. ;;; Generate extensions for radicals of polynomials
  175. ;;; Currently this does not split previously defined radicands.
  176. ;;; It will as soon as expression rework is added.
  177. (define (make-radical-ext p r)
  178.   (set! p (licit->polxpr p))
  179.   (let ((prest #f)
  180.     (pegcd #f)
  181.     (radrest #f)
  182.     (en #f)
  183.     (e (member-if (lambda (e) (equal? p (cadr e))) radical-defs)))
  184.     (cond (e (if (divides? r (length (cddr (car e))))
  185.          (radpow (car e) r)
  186.          (var->expl (make-rad-var p r))))
  187.       ((eqv? 0 p) (var->expl (make-rad-var p r)))
  188.       ((begin (set! e (member-if (lambda (rule)
  189.                        (set! en (cadr rule))
  190.                        (set! pegcd (poly:gcd en p))
  191.                        (not (eqv? 1 pegcd)))
  192.                      radical-defs))
  193.           e)
  194.        (set! prest (poly:/ p pegcd))
  195.        (set! radrest (poly:/ en pegcd))
  196.        (if (and (eqv? 1 radrest) (divides? r (length (cddr (car e)))))
  197.            (app* $1*$2 (make-radical-ext prest r) (radpow (car e) r))
  198.            (var->expl (make-rad-var p r))))
  199.       (else (var->expl (make-rad-var p r))))))
  200.  
  201. (define (radpow radrule r)
  202.   (univ:monomial 1 (quotient (length (cddr radrule)) r) (car radrule)))
  203.  
  204. ;;;    Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
  205.